Fast and accurate computation of Jacobi expansion coefficients of analytic functions

The computation of spectral expansion coefficients is an important aspect in the implementation of spectral methods. In this paper, we explore two strategies for computing the coefficients of polynomial expansions of analytic functions, including Chebyshev, Legendre, ultraspherical and Jacobi coefficients, in the complex plane. The first strategy maximizes computational efficiency and results in an FFT-based O(N logN) algorithm for computing the first N spectral expansion coefficients. This strategy is stable with respect to absolute errors and recovers some recent algorithms for the computation of Legendre and ultraspherical spectral coefficients as special cases. The second strategy maximizes computational accuracy. We show that an optimal contour in the complex plane exists for each Chebyshev expansion coefficient. With these contours Chebyshev coefficients can be computed with small relative error, rather than the usual small absolute error, at a cost of increasing the computational complexity. We show that high accuracy is maintained even after repeated differentiation of the expansion, such that very high order derivatives of analytic functions can be computed to near machine precision accuracy from their Chebyshev expansions using standard floating point arithmetic. This result is similar to a result recently obtained by Bornemann for the computation of high order derivatives by Cauchy integrals. We extend this strategy to the accurate computation of Jacobi coefficients.